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Abstract 

This article presents a comparative study of three different types of estima- 
tors used for supervised linear unmixing of two MEx/OMEGA hyperspectral 
cubes. The algorithms take into account the constraints of the abundance 
fractions, in order to get physically interpretable results. Abundance maps 
show that t he Bayesian maximum a post e riori p robability (MAP) estimator 



proposed in iThemelis and Rontogiannis 



(120081 ) outperforms the other two 



schemes, offering a compromise between complexity and estimation perfor- 
mance. Thus, the MAP estimator is a candidate algorithm to perform ice 
and minerals detection on large hyperspectral datasets. 
Keywords: Hyperspectral imagery, supervised unmixing, OMEGA data, 
Mars Express 
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1. Introduction 



The surface of Mars is currently being imaged with a combination of 
high spectral and spatial resolution. This gives the ability to detect and 
map chemical components on the Martian surface and atmosphere more ac- 
curately than before. Spectr al unmixing (SU) is one of t h e techniques cur 



(|2008h : 



Schmidt et al. 



rently used for this pu r pose, Keshava and Mustardl (120021 ); 



Moussaoui et al 



(120101 ). SU is the procedure by which the measured 
spectrum of a mixed pixel is decomposed into a number of constituent spec- 
tra, called endmembers, and the corresponding fractions, or abundances, 
that indicate the proportion of the presence of each endmember in the pixel, 
Keshava and Mustardl (120021 ). Linear SU, which adopts the hypothesis that 
the spectrum of a mixed pixel is a linear combination of its endmembers' 
spectra, is more commonly used in practice. Based on a physical interpre- 
tation, two hard constraints are imposed on the abundance fractions of the 
materials in a pixel; they should be non-negative and sum to one. 

Several SU techniques for the unmixing of OMEGA (Observat oire pour la 
Miner alogie, 1' Eau, les Glaces et 1' Activite hyperspectral images) . iBibring et al 
( 120041 ). have been recently proposed in the bibliography. These techniques 
can be categorized into unsupervised, where a special procedure is first ex- 
ecuted to get the endmembers' spectral signatures from the image, and su- 
pervised, where a priori knowledge of the image endmembers is available. 
A recent example of an unsu pervised technique is the Bayesian source sepa- 



ration method, developed in 



Schmidt et al. 



(120101 ) . This technique is based 



on a Gibbs sampling scheme to perform Bayesian inference. Due to its high 
computationally complexity, a special implementation strategy is developed 
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for its application to a complete OMEGA image data set, iSchmidt et al. 
(120101 ). B and ratio is th e most commonly used supe r vised technique to dete ct 



minerals (IBibring et al 



(120041 ): 



Bibring et al. 



( 120041 ): 



Langevin et al 



(12005( 1 ). 



However some multiple endmember 
been proposed such as MELSUM, 



inear spectra l unmi xing algorithms have 



Combe et al 



(120081 ). MELSUM uses a 



reference library containing various spectral signatures of minerals (used as 
endmembers), and is based o n the classical spectral mixture analysis (SMA) 



algorithm. In 



Kanner et al. 



( 120071 ). a modified Gaussian model (MGM) has 
been exploited to estimate the fractional abundances of a compositionally 
diverse suite of pyroxene spectra in the martian surface. A wavelet based 
method has also been appli e d for the unm i xing o f hyperspectral martian 



data in 



Gendrin et al. 



(120061 ): 



Schmidt et al. 



(120071 ). 



In this paper, we focus on the problem of supervised SU. Our main ob- 
jective is to estimate the abundances of the endmembers that are present 
in two OMEGA images, subject to the non-negativity and sum-to-one con- 
straints. In the following, three different supe rvised unm i xing a lgorithms are 



considered, namely the ENVI- SVD method, 



Boardmanl (119891 ). a quadratic 



programming (QP) technique, IColeman and Lil (119961 ). and a recently pro- 
posed Ba yesian maximum a posterior i prob ability soft-constraint (MAPs) 



estimator, 



Themelis and Rontogiannid (120081 ) . These algorithms are applied 



on two different hyperspectral OMEGA data sets, and they are evaluated 
through their corresponding abundance maps. The endmember reference 
spectra used for supervised unmixing are either extracted from the image 
itself, or selected from a spectral library of pure minerals. The experimen- 
tal results show that the MAPs estimator results in abundance values that 
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satisfy the constraints of the problem and provides a compromise between 
the performance of the QP technique and the complexity of the ENVI-SVD 



method. An earlier version o 
Planetary Science Congress, 



this paper was presen ted at the 2010 European 



Themelis et al. 



(boid ) 



2. Linear spectral unmixing techniques 

Before we present the unmixing techniques considered in this paper, 
a short description is presen ted here on the linear mixing model (LMM), 
Keshava and Mustard! (120021 ). which is assumed in all three of them. In a 
hyperspectral image, each pixel is represented by a L- dimensional vector 
y, where L is the number of the available spectral bands. The elements of 
y correspond to the reflectance measured at the respective spectral bands. 
The LMM assumes that the received pixel's spectrum is generated by a linear 
combination of endmembers' spectra. Suppose that the spectral signatures 
of p materials that may exist in the image are available. Then, y can be 
expressed by the following linear regression model: 



y = <frx + n, 



(1) 



where * = [(j)^ . . . (j> J G M+ Xp , is the mixing matrix containing the end- 
members' spectra (L-dimensional vectors z = 1, 2, . . . ,p), x is a p x 1 
vector with the corresponding abundance fractions, and n is a L x 1 additive 
noise vector. 

Adopting the linear model in (JTJ), three different unmixing algorithms 
are applied to the O MEGA data sets: i) a singular value decomposition 
method (ENVI-SVD), iBoardmanl (119891). available in the E NVI image pro- 



cessing software ii) a QP technique, 



Coleman and Li 



(119961 ). available in the 



Matlab environment, and hi ) a rec ently proposed Bayesian MAPs estimator, 
Themelis and Rontogiannisl ( 120081 ). These algorithms are briefly described 
in the following subsections. 



2.1. ENVI-SVD 

ENVI-SVD is a constrained least squares approach to the unmixing prob- 
lem. Using the singular value decomposition (SVD) algorithm, the pseudo- 
inverse of the mixing matrix <fr is computed. Then, the abundance fractions 
are easily estimated by multiplying the pseudo-inverse matrix with each im- 
age pixel's spectral vector. The advantage of this method is its low compu- 
tational complexity, since the pseudo-inverse matrix is computed only once 
as a preprocessing step, and is then applied to all image pixel vectors. As far 
as the sum-to-one constraint is concerned, it is imposed to the problem using 
an extra (weighted) equation to the linear system of equations ([1]). However, 
ENVI-SVD does not take into account the non-negativity of the abundances, 
which can result in negative abundance values that have no physical meaning. 

2.2. Quadratic programming technique 

This quadratic programming technique is a reflective Newton method, 
which minimizes the quadratic function of the least squares error of the 
unimixing problem, subject to the sum-to-one and non-negativity constraints 
imposed on the abundances. The QP technique is based on an iterative op- 
timization scheme, which has to be repeated separately for each image pixel 
vector. Keeping in mind that a hyperspectral image may be composed of 
thousands of pixels, solving a separate optimization problem for each pixel 
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adds up to the computational complexity of the method. An implementa- 
tion of the quadratic programming technique is available in the optimization 
toolbox of Matlab. 



2.3. MAPs estimator 



The MAPs estimator proposed in iThemelis and Rontogiannisl (120081 ) is 



a Bayesian estimator specifically designed to address the inverse problem 
of supervised hyperspectral unmixing. In a Bayesian framework, appropri- 
ate prior distributions are assigned to the unknown parameters of the esti- 
mation problem, which usually reflect the paramete rs' natural characteris- 



tics. The Bayesian approach of 



Benavoli et al. 



( 120071 ) is adopted in which a 



Gaussian distribution is used as a prior for the abundance vector x and the 
MAP estimator is then utilized. By exploiting the symmetry of the prob- 
lem's convex constraints, the parameters of the Gaussian posterior distribu- 
tion (i. e., the mean and the covar i ance m atrix) can be expressed in closed 



forms, 



Themelis and Rontogiannisl ( 120081 ). Due to the statistical nature of 



the Bayesian estimator, the constraints are not explicitly imposed to the esti 
mated parameters. To alleviate this, the final step of the algorithm is a pro 
jection of the MAP estimation point on the polyt ope of constraints, providing 



its nearest estimate that satisfies the constraints. 



Themelis and Rontogiannis 



(120081 ). This algorithm has substantially lower complexity than the QP tech- 



nique, since it relies on the computation of simple closed-form expressions. 
3. Discussion 

The previously described algorithms are applied to two different OMEGA 
data cubes: (a) a scene of Mars' South Polar Cap, and (b) a Syrtis Major 
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observation. The OMEGA instrument is a spectrometer on board ESA's 
Mars Express satellite, which provides hyperspectral images of the Mars sur- 
face, with a spatial resolution from 300m to 4km, 96 wavelength channels 
in the visible band and 256 wavelength channels in the near infrared band, 
( 12004] ) . OMEGA uses three different detectors, with spectral 



Bibring et al 



resolutions about 7.5nm in the 0.35 — 1.05/im wavelength range (visible and 
near infrared channel or VNIR), 14nm between 0.94 and 2.70/iin (short wave 
infrared channel or SWIR) and an average of 21nm from 2.65 to 5.2/iin (long 
wave infrared channel or LWIR), respectively. The two hyperspectral data 
sets, the reference spectra used for each image and the SU results obtained 
from the application of the three methods are analytically described in the 
following sections. 

3.1. South Polar Cap image cube 

This data set consists of a single hyperspectral data cube obtained by 
looking towards the South Polar Cap of Mars in the local summer (Jan. 
2004). The data cube is made up of two channels: 128 spectral planes from 
0.93 to 2.73/iin with a resolution of 14nm and 128 spectral planes from 2.55 
to 5.11/xm with a resolution of 21nm. Noisy bands were excluded, and 156 
out of the 250 initial bands were finally utilized in the region from 0.93 
to 2.98/im to avoid the thermal emission spectral range. The linear model 
mixing matrix consists of the following three reference spectra: (a) CO2 ice 
(synthetic data with grain size = 100mm), (b) H 2 ice (synthetic data with 
grain size = 10/rni), and (c) dust, whi c h wer e all detected a priori using 



the Wavanglet method of 



Schmidt et al 



cussed in the first OMEGA publication, 



12007). These endmembers are dis- 



Bibring et al 



(j2004h . and are also 



verified bv lSchmidt et al.l (12 010) 



tion method of 



Moussaoui et al. 



using the Bayesian positive source separa- 



(120081 ) . The respective spectral signatures 
of the endmembers are shown in Fig. [TJ 

The abundance maps for each endmember resulting after the application 
of the three estimators are displayed in Figs. [2]-|U As shown in these figures, 
the ENVI-SVD abundances do not satisfy the non-negativity constraint, e.g., 
regarding the CO2 endmember, the minimum computed abundance value 
is —9.9 x 10~ 2 . Notice also that the abundance values calculated by QP 
and MAPs are in full agreement, they share the same scale and are quite 
different from those obtained by ENVI-SVD. This shows that the MAPs 
estimates provide reliable information about the abundances. It also adds 
up to the fact that the MAPs estimator has almost similar performance with 
the QP algorithm in simulation s cenarios with synthetic data, as shown in 
Themelis and Rontogiannisl (120081 ) . 



3.2. Syrtis Major image cube 

This data set consists of a single hyperspectral data cube of the Syrtis Ma- 



jor region, whic h contains wel 
mafic minerals 



Mustard et al. 



- identi fied areas with very strong signatures of 
( 120051 ). The data cube consists of 109 spectral 
bands out of the 128 original wavelengths of the SWIR detector. The spatial 
dimensions of the cube are 366 x 128 pixels. The OMEGA observations have 
been calibrated for known instrument artifacts and for atmospheric CO2. 
The cube has been radiometrically corrected using the standard correction 
pipeline (SOFT06) and the atmospheric gas tra nsmission has b e en em pir- 



ically corrected using the volcano scan method, 



Langevin et al. 



(120051 ). It 



is well known that OMEGA can identify pyroxene and olivine; it discrimi- 



nates between the high-calcium pyroxenes (HCPs, e.g ., clinopyroxen e s) an d 



Bibring et al. 



()2005|) 



low-calcium pyroxenes (LCPs, e.g., orthopyroxenes) 

In the scene under investigation, we utilized three e ndmembers whi c h hay e 

J2OO5L 



previously been identified to be present in the image, 



Mustard et al. 



namely, (a) Hypersthene, (b) Diopside (c) Fayalite. These are all l abora 



tory reference spectra, which have also been used in iSchmidt et al.l (120111 ) 
for supervised unmixing. It is interesting to note that the last two endmem- 
bers have been retrieved using CRISM multispectral observations. Their 
respective spectral signatures are displayed in Fig. |5j Three more arti- 



fact endmember spectra are utilized, specifically two neutral spectra 



ponents (flat lines at 10 4 and 1), and a slope line, as in 



Schmidt et al. 



com- 



Le Mouelic et al 



( 120091 ) 



(|201l|); 



The abundance maps obtained from the application of the three methods 
to the Syrtis Major hyperspectral scene are shown in Figs. |6]-[8j Each Figure 
illustrates the corresponding abundance map of a single endmember, as it is 
estimated by all three methods. As a reference, the abundance maps of all 
three endmembers using the band ratio method are also shown in Fig. |9j 
The band ratio and bands depth estimation methods are commonly used to 

(120051 ) . This method is valid 



detect minerals on OMEGA data 



Bibring et al 



only if at least two wavelength channels can be identi fied as affected by on 



one particular mineral. Following the methodology of 
we used four band ratios: 



Schmidt et al. 



(2011), 



y 



Index(Olivine) = b2.39/bl.06 



(2) 
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Index(opx) = 1 - bl.84/((1.84 - 1.25) * bl.25 + (2.47 - 1.84)* 



b2.47) * (2.47 - 1.25) 



(3) 



Index(cpx) = 1 - bl.85/((2.32 - 1.85) * b2.32 + (2.56 - 2.32) 



*b2.56) * (2.56 - 1.85)) 



(4) 



where the wavelength band "bl.84" stands for the band at 1.84 microns, 
Orthopyroxenes (Hyperstene) are noted as "opx" and Clinopyroxenes (Diop- 
side) are noted as "cpx". The differences between the three methods are 
again prominent. Both the MAPs and the QP methods provide consistent 
results, as far as the endmembers Fayalite and Hypersthene is concerned. 
In addition, for both Fayalite and Hypersthene, ENVI-SVD returns negative 
abundance values and its resulting maps substantially deviate from the other 
two methods. As far as Diopside is concerned, the MAPs estimator seems to 
produce a slightly different abundance map in comparison to the other two 
methods. However, it can be easily verified by comparing Figs. [6] - [8] and 
[9] that the abundance maps of the MAPs estimator are in better agreement 
with those obtained using the band ratio method, compared to the other two 
methods. Thus, it can be argued that MAPs provides more reliable results 
than ENVI-SVD and QP. 

4. Conclusions 

In this paper, we have presented a comparison of three different super- 
vised spectral unmixing methods (Bayesian MAPs estimator, ENVI-SVD, 
QP), on the basis of two different OMEGA hyperspectral data sets. As op- 
posed to iterative algorithms or Markov Chain Monte Carlo methods, e.g. 
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Schmidt et al.l (120101 ) , commonly used for the constrained inverse problem of 
abundance estimation, the computational complexity of the MAPs estimator 
is much lower. Specifically, for the Syrtis Major dataset, the running time 
of the MAPs algorithm was 4.7 sees in a roughly optimized Matlab imple- 
mentation, while the QP needed 26.4 sees (both algorithms were run on a 
2.4-Ghz Intel Core 2 CPU). In addition, as verified by experimental results, 
the performance of the two methods is approximately equal. Therefore, the 
MAPs estimator seems to offer the best compromise between estimation per- 
formance and complexity among the three algorithms, and it is thus a serious 
candidate for spectral unmixing of hyperspectral data in Planetary Sciences. 
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Figure 1: Reference spectra of the South Polar Cap OMEGA image. The 
available endmembers are: (a) OMEGA typical dust materials with atmo- 
sphere absorption, (b) synthetic C02 ice with grain size of 100 mm, (c) 
synthetic H20 ice with grain size of 100 microns. 
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Figure 5: Reference spectra of the Syrtis Major OMEGA image. 
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Figure 9: Abundance maps of the endmembers (a) Hypersthene, (b) Diop- 
side, and (c) Fayalite in the Syrtis Major scene. The abundances are esti- 
mated using the band ratio method. 
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